#!/usr/bin/perl
# 2013/10/16 ver1. extract feature and make BED file

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 8 or $ARGV[0] eq '--help'){
	die "Usage : $0 -a [contig for chr1] -b [contig for chr2] -c [contif for chr3] -d [output directory]\n";
}

my %opt;
getopts("a:b:c:d:", \%opt);
my $FILE_chr1 = $opt{a};
my $FILE_chr2 = $opt{b};
my $FILE_chr3 = $opt{c};
my $DIR = $opt{d};

# output file
my $FILE_tRNA = $DIR . '/tRNA.bed';
my $FILE_LTR = $DIR . '/LTR.bed';
my $FILE_5SrRNA = $DIR . '/5SrRNA.bed';
my $FILE_CDS = $DIR . '/CDS.bed';
my $FILE_ribosome = $DIR . '/ribosomal_gene.bed';

my $fh_tRNA = IO::File->new($FILE_tRNA, 'w') or die "cannot write $FILE_tRNA: $!";
my $fh_LTR = IO::File->new($FILE_LTR, 'w') or die "cannot write $FILE_LTR: $!";
my $fh_5SrRNA = IO::File->new($FILE_5SrRNA, 'w') or die "cannot write $FILE_5SrRNA: $!";
my $fh_CDS = IO::File->new($FILE_CDS, 'w') or die "cannot write $FILE_CDS: $!";
my $fh_ribosomal = IO::File->new($FILE_ribosome, 'w') or die "cannot write $FILE_ribosome: $!";

&parseFile($FILE_chr1, 'I');
&parseFile($FILE_chr2, 'II');
&parseFile($FILE_chr3, 'III');

$fh_tRNA->close();
$fh_LTR->close();
$fh_5SrRNA->close();
$fh_CDS->close();
$fh_ribosomal->close();

sub parseFile{
	my ($FILE, $chr) = @_;
	my $TYPE = '';
	my $block = 0;
	my @size;
	my @starts;
	my @ends;
	my @starts_conv;
	my $direction;
	my $gene = '';
	my $systematicID = '';
	my $product = '';

	my $fh_in = IO::File->new($FILE) or die "cannot open $FILE: $!";
	while($_ = $fh_in->getline()){
		s/\r?\n//;
		if(/^SQ/){
			if($gene eq ''){
				$gene = $systematicID;
			}
			if($TYPE eq 'tRNA'){
				$fh_tRNA->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($TYPE eq 'rRNA'){
				$fh_5SrRNA->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($TYPE eq 'LTR'){
				$fh_LTR->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], "LTR", $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($product eq 'Ribosomal gene'){
				$fh_ribosomal->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($TYPE eq 'CDS'){
				$fh_CDS->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
		}
		unless(/^FT/){
			next;
		}

		my $body = substr($_, 21);
		my $tmp_TYPE = substr($_, 5, 15);
		$tmp_TYPE =~ s/\s//g;

		if($tmp_TYPE ne ''){
			# check and output previous data
			if($gene eq ''){
				$gene = $systematicID;
			}
			if($TYPE eq 'tRNA'){
				$fh_tRNA->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($TYPE eq 'rRNA'){
				$fh_5SrRNA->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($TYPE eq 'LTR'){
				$fh_LTR->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], "LTR", $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($product eq 'Ribosomal gene'){
				$fh_ribosomal->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}
			if($TYPE eq 'CDS'){
				$fh_CDS->printf("%s\t%d\t%d\t%s\t0\t%s\t%d\t%d\t0\t%d\t%s\t%s\n",
				$chr, $starts[0], $ends[-1], $gene, $direction,  $starts[0], $ends[-1], $block, join(",", @size), join(",", @starts_conv));
			}

			# initialize
			$TYPE = $tmp_TYPE;
			$block = 0;
			@size = ();
			@starts = ();
			@starts_conv = ();
			@ends = ();
			$gene = '';
			$systematicID = '';
			$product = '';
			$direction = '+';
			if(/complement/){
				$direction = '-';
			}

			# check location
			while(/(\d+)\.\.\D?(\d+)/g){
				my $s = $1;
				my $e = $2;
				push @starts, $s;
				push @ends, $e;
				push @size, ($e - $s);
				push @starts_conv, ($s - $starts[0]);
				$block++;
			}
		}
		if($body =~ m|/gene="([\w\.]+)"|){
			my $tmp_gene = $1;
			if($gene eq '' or $gene =~ /\./){
				$gene = $tmp_gene;
			}
		}
		if($body =~ m|/product=".+ribosomal protein|){
			if($body !~ /mitochondrial/){
				$product = 'Ribosomal gene';
			}
		}
		if($body =~ m|/systematic_id="([\w\.]+)"|){
			$systematicID = $1;
		}
	}
	$fh_in->close();
}






